Exact calculation of the energy contributions to the T = random-field Ising model 
with metastable dynamics on the Bethe lattice 
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We analyze the energy terms corresponding to the spin-spin exchange ^ SiSj and spin-random 
field coupling ^ Sihi of the zero temperature random-field Ising model on the Bethe lattice driven 
by an external field with metastable dynamics. Exact results are calculated as a function of the 
standard deviation of the disorder a and the coordination number z, and compared with numerical 
simulations on random graphs for z = 4, for which a disorder-induced transition takes place. 
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I. INTRODUCTION 

Hysteresis and metastability are intriguing phenom- 
ena with implications in both fundamental and applied 
physics 0. They arise as a consequence of the exis- 
tence of internal energy barriers that cannot be overcome 
by the system. A particularly interesting case is the so 
called "rate-independent" hysteresis, for which metasta- 
bility is not a consequence of the fast driving rate but 
of the "athermal" character of the system. The energy 
barriers are so high compared with thermal fluctuations 
that effectively the system behaves at zero temperature, 
following a reproducible and deterministic metastastable 
path when the external field is varied. Many experimen- 
tal situations are known to be well approximated by this 
extreme case. 

Several models have been useful for the characteriza- 
tion and description of rate independent hysteresis in 
different experimental systems. A very interesting mi- 
croscopic model is the Random Field Ising model at 
T = with metastable dynamics (T = RFIM) 0. 
It contains the three essential competing ingredients for 
the occurence of "athermal" hysteresis: a ferromag- 
netic nearest-neighbour (n.n.) interaction term favour- 
ing long range order, a local energy term associated with 
quenched disorder, and a term describing the coupling of 
the magnetization to an external driving field. Although 
the model is formulated using a magnetic language it can 
be translated easily to other systems displaying athermal 
hysteresis. 

During the last decade the T = RFIM model has 
been used as a prototype for the understanding of many 
properties associated to hysteresis: return point mem- 
ory, congruency, distribution of metastable states 0, Q , 
demagnetization process Q , disorder induced critical 
points 8, 9] and power-law distribution of the magne- 
tization avalanches (Barkhauscn noise) 0, ^| . Most of 



these results are based on numerical simulations on finite 
lattices or on mean field analysis. Interestingly, however, 
non trivial analytical solutions can be obtained for the 
particular case of the T — Q RFIM on Bethe lattices 
with coordination number z. The main hysteresis loop 
was solved seven years ago by Shukla, Dhar and Sethna 
[r3.[T3l |. More recently partial loops fl^,'!^ and trajecto- 
ries starting from states with "quenched" spins 16] have 
also been deduced . Here we present the explicit com- 
putation of the three energetic terms of the Hamiltonian 
TL for a generic values of a and z 19]. Analytical results 
for all the contributions to TL give insight on the singular 
behaviour of the system at Uc, and allows for a deeper 
understanding of the energy balances in the hysteresis 
loop. 

The paper is organized as follows. In section ^ we 
summarize the details of the model. In section IIIII we 
solve the different terms in the Hamiltonian. In section 
llVl we present the results for the z = 4 case and compare 
with numerical simulations on random graphs. Finally, 
in section we study the energy dissipation. 



II. MODEL 

The RFIM is defined on a Bethe lattice with N sites 
and coordination number z. On each lattice site we define 
a spin variable Si which takes values ±1. The Hamilto- 
nian (magnetic enthalpy) reads: 



-Ue + Ud- HM 



where H is the external driving field, 



N 



(1) 



(2) 



is the magnetitzation, 
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U,{{Si])^-jY,S,S, (3) 
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is the ferromagnetic exchange energy extending over all 
n.n. pairs, and 



N 



(4) 



of the random fields of the three terms in the Hamiltonian 
gives: 



-HM) 



-Hm 



(7) 



accounts for the energy interaction with quenched dis- 
order. The random fields {hi} are independent and dis- 
tributed according to a Gaussian probability density cen- 
tered around zero: 



fih^) = 



1 



2tt(t 



(5) 



where a is the standard deviation of the random fields 
and controls the amount of disorder in the system. 

For the analysis of metastability and hysteresis loops 
we use a 1-spin-flip local relaxation dynamics. This 
is the standard choice used in previous studies of the 
metastable T = RFIM 0]: each spin Si flips individ- 
ually according to the sign of its local field Fi given by: 



H + h, 



(6) 



where the first sum extends over the z neighbours of Si. 
The complete lower branch of the hysteresis loop is ob- 
tained by adiabatically increasing H from —oo (RI ~ 
-N) to +00 (M = -l-A^). 

In order to check the analytical results that will be 
presented we have also performed numerical simulations 
on random graphs with coordination number z = 4 and 
a range of sizes from N — 10"* to = 10^ . It is known 
that in the thermodynamic limit such numerical simula- 
tions agree with the analytical results for Bethc lattices 
[l3| . For the simulations we start with a value of H neg- 
ative enough so that the unique stable state is given by 
all the spins Si = —1. Wc increase the external field H 
until Fi vanishes on a certain spin. The spin is then re- 
versed keeping H constant. This reversal may destabilize 
some of the neighbouring spins, which are then reversed 
simultaneously (parallel updating). This is the begin- 
ing of an avalanche. The avalanches proceed until a new 
metastable situation with all the spins Si aligned with 
their respective local fields Fi is reached. We can then 
continue increasing the external field H. 

In the figures presented below numerical simulations 
correspond always to averages over several (~ 10) real- 
izations of the random graph and many realizations of 
the random fields 1000). 



III. GENERAL SOLUTION FOR 
COORDINATION NUMBER z 

Our goal is the computation of the different energetic 
terms in the Hamiltonian. The average over realizations 



{Ue) _ ^ j,^^ 



N 



{hiSi 



(8) 



(9) 



The averages on the right hand sides can be computed 
as follows: 

(5.) = ^5,P(5.) = mi)-^(-l) (10) 

where P{Si) is the probability for a spin to take a value 
±1, 

{S,S,) = J2 S^,S,PiS,,S,) = 
{s.,s,} 

= P(+1,+1) + P(-1,-1)-2P(+1,-1) 

(11) 

where P{Si, Sj) is the probability for a nearest neighbour 
pair to be in the state {Si,Sj). This probability satisfies 
P(+l,-l) = P(-l,-|-l). Finally 

/ + 00 
dhiP{h^,Si)hiSt = 

/+00 r+oo 
dh,P{h,,+l)h,- / dh,P{h,,^l)h, (12) 
-oo J —oo 



where P{hi, Si) is the probability for a site i of having a 
random field within (hi, hi + dhi) and a spin with state 
S,. 

At this point we follow Ref. 0. We define the proba- 
bility P{Si\n) for a spin being in state Si given a certain 
environment of nearest neighbours. This environment is 
fully characterized by the variable n (0 < n < z) which 
accounts for the number of neighbours in state -|-1 (see 
Fig.UIa)). Clearly, 

/+ca 
dhJih,) = 

1 , ( -J(2n- z) - H 
= -erfc < ^ 

2 I V2a2 

and 

P(5, = -l|n) = l-P(S, = -fl|n) 
From Bayes formula one can write: 



p(5,) = E^(^^)^(^»i") 



(13) 
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where P{n) is the probabihty for a site having an envi- 
roment with n neighbouring spins in the state +1 (see 
Fig. n^a)). According to Ref. [13, 

PH= (^^)p*"(l-P*)(-") (14) 

where P* results from defining the conditional probabil- 
ity that a spin is +1 given that its parent spin (according 
to the hierarchy of the Bethe lattice) is down, and its de- 
scendent spins are relaxed. In a site deep enough inside 
a very big lattice (thermodynamic limit), this probabil- 
ity tends to be homogeneus and P* is given by the self 
consistent equation 

2-1 

P* ^™(^ ~ P*y-^-''P{S^ = +l|n) (15) 

Note that this equation implicitely contains the infor- 
mation on the fact that we are increasing the field 
monotonously from the negative fully saturated state. 

By numerically obtaining P* from H15|) and using H14|) 
and H13|) . the averaged magnetization in (|10|l can be 
obtained. This allows to calculate the averaged lower 
branch of the hysteresis loop. 

In order to compute the terms 1)11(1 and (|12|l we apply 
Bayes formula again and write 

P(S„5,) =^^P(/,r)P(^„5,|;,r) (16) 

/=1 r=l 



P(/i„^,) = I]^H^(/i.,^.|'^) , (17) 

n=l 

where P{Si, Sj\l^r) and P(hi, Si\n) are conditional prob- 
abilities given a certain environment and P{l,r) is the 
probability for a pair having an environment with / spins 
in the state +1 in the left neighbourhood and r spins in 
the state +1 in the right neighbourhood (see Fig. QJb)). 
This is the generalitation of P(n) for the description of 
the environment of a pair of spins. 

A. Calculation of P{Si, Sj) 

The calculation of P{Si,Sj) starts by generalizing 
Eq. (jUI) to 

P(;,r) = 5(Z,r)P*'(l -P*)(^-1-')P*''(1 -P*)(^-i-'-) ^ 
= 5(Z,r)P*'+''(l-P*)(2--2-r-;) (18) 

where 

The conditional probabilites for a pair of spins given 
a certain environment, P{Si, Sj\l^r), can be written as 




I spins 4-1 r spins +1 

FIG. 1: Schematic representation of the enviroment of (a) a 
single site and (b) a pair. The variables n, I, r account for 
the number of spins -1-1. 
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FIG. 2: The {hi, hj) plane divided in different areas (different 
color) corresponding to the final state of spins {Si,Sj) for a 
fixed enviroment Z = l,r = 1. These areas correspond to the 
domains D of integration in Eq. 1201 . 



double integrals of the random field distribution on a 
certain domain in the hi — hj plane, i.e. 

P{S,,Sj\l,r) = J J^f{h,)f{hj)dh,dh, (20) 

Figure |2] shows an example of such domains correspond- 
ing to the environment I — 1 and r — I. The plots corre- 
sponding to other environments are obtained by transla- 
tion of this case. Therefore: 
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P{-l,-l\l,r) = [1 - P(5, = +1|0] [1 - P{S^ = +l|r)] 

(21) 



r) = P{S^ = [1 - = +l|r + 1)] 

(22) 

P(-l, r) = [1 - P(5, = +1|; + 1)] P(5, = +l|r) 

(23) 

P(+l,+l|/,r) = P(5, = +l|/ + l)P(5, = +l|r) + 
+P(5. = +1|0 [P{S^ = +l\r + 1) - P(5, = +l|r)](24) 

From these equations and using p8|l and one can 
obtain the correlation (SiSj) needed for the computation 
of Eq. (|TT]l. 

B. Calculacion of P{hi,Si) 

The computation of P{hi, Si) is straightforward noting 
that P{hi, Si\n) can be computed as tails of the gaussian 
distribution f{hi) as: 

Thus, using p4|) and p7|) one can compute {hiSi) needed 
for the computation of Eq. (fT^ . 

IV. NUMERICAL SOLUTION FOR THE CASE 

The case with coordination number z = 4 is inter- 
esting because it is known to display a disorder in- 
duced phase transition between smooth hysteresis loops 
for (T > ac — 1.78125 and discontinuous hysteresis loops 
for (T < (Tc. We present the results obtained by nu- 
merically solving the equation of the previous section 
for 2 = 4. In particular the real roots of Eq. (|15|l are 
found with a bisection method restricted to the interval 
< P* < 1. 

Figure|21shows the m-H diagram corresponding to four 
different amounts of disorder a. Note that the data rep- 
resented corresponds only to the lower branch of the hys- 
teresis loop for increasing the external field H. 




H H 

FIG. 3: Lower branches of the hysteresis loop corresponding 
to z = 4 and different amounts of disorder. Exact results 
(continuous lines) are compared with numerical simulations 
(dotted lines). Simulations correspond to a random graph 
with N = 10^ and averages over 1000 different realizations of 
disorder 

observed in Fig. O for a < (Tc- Nevertheless, numerical 
simulations show that only one of the roots has physical 
meaning. In the case of increasing field, only the lower m 
branch is obtained in the simulations, and the disconti- 
nuity occurs at H2 where this lower branch of the s-shape 
curve joints the intermediate branch and disappears. To 
our knowledge there is no clear physical explanation for 
this fact and, a priori, the jump could occur at any field 
in the range Hi-H2. A stability criterium would be de- 
sirable. In this respect we note that in a very recent 
paper it has been speculated that the s-shape curve 
is related to the boundary of the density of 1-spin-flip 
metastable states 17]. It is also important to notice that 
the result of numerical simulations depends on system 
size, as shown in Fig.^ Only in the thermodynamic limit 
[N — > 00) the data from numerical simulations would fol- 
low exactly the lower branch of the theoretical curve up 

to H2. 



Figure |31 shows the behaviour of Ue corresponding to 
the same four cases as in Fig. O Below CTc the three 
roots of Eq. (|15|) generate a lace function. Numerical 
simulations follow continuously one of the roots until H2 
where they jump to the lower exchange energy branch. 
Note that there is an intermediate field (crossing point 
of the lace) where two of the roots correspond to the same 
value of Ue- 



The disconinuity in the m{H) branch, as was pointed 
out in Ref. > arises from the fact that the solution of 
equation IjlSfl is trivalued for cr < CTc in a certain field 
range Hi < H < H2- The three roots of Eq. H15|l gener- 
ate the s-shape curve in the m — H diagram than can be 



Figure shows the behaviour of Ud [l^- A similar lace 
shape is obtained. We note, that the crossing points 
(for a < ac) are different from in the Ue curve. As a 
tends to ac the fields Hu H2, H3, and all tend to the 
critical value He = 1 |13| . 
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FIG. 4: Examples of finite size dependence of the lower 
branch of the hysteresis loop. Symbols correspond to sim- 
ulations on random graphs with increasing sizes as indicated 
by the legend. The continuous line corresponds to the exact 
solution. 
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FIG. 5: Exchange energy behaviour corresponding to the 
same cases as in Fig. |3| Lines correspond to numerical solu- 
tion of the exact equations and dots correspond to numerical 
simulation on random graphs. 
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FIG. 6: Disorder-coupling energy corresponding to the same 
cases as in Fig. |3 Lines correspond to numerical solution of 
the exact equations and dots correspond to numerical simu- 
lation on random graphs. 




FIG. 7: Behaviour of the Hamiltonian 7i as a function of 
the external field H corresponding to the same four cases as 
the previous figures. Numerical solution of the exact equa- 
tions are shown with a continuous line, whereas numerical 
simulations of random graphs are shown with a dotted line. 



In Fig. [71 we show the behaviour of the total Hamilto- 
nian (magnetic enthalpy) TL, corresponding to the same 
four cases as before. The plots show that, in the trivalued 
region, the numerical simulations choose the branch with 
maximum 7i, which is clearly a non-equilibrium path. 



V. ENERGY DISSIPATION 

The fact of having obtained separately m, U^, and U^^ 
as functions of H, allows an exact computation of the 
energy dissipated by the system Q. It can be calculated 



by integration along the non-equilibrium path [ig: 

Q = AU- J HdM, (27) 

where AU — AUe + AUd (internal energy difference 
between the initial and final states). Starting from 
H — — oo, the path integral over the trajectory com- 
puted in the previous section can be written in the form: 

Q(-oo H') U{H') 1 ^ , 

— - + -zJ- H dm. (28) 

Figure IHl shows the result obtained from this expression 
for z = 4 and different values of a. For those cases for 
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which a < ac the energy dissipated has been computed 
for three different trajectories (indicated by dashed or 
continuous hues), all of them compatible with the theo- 
retical (numerical) solution of m{H), Ue{H), and Ud{H): 
the first one assumes that the transition to the upper 
branch of m{H) takes place a.t H — Hi, the second one 
at an intermediate value of the field H = {Hi + H2)/2, 
and the third one aX, H — H2- Of these trajectories, we 
verify that the one jumping at H2 (the trajectory chosen 
by the system in the numerical simulations) is the one of 
maximum energy dissipation [Q most negative). 
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FIG. 8: Energy dissipated along the metastable path -co — > 
H, computed from Eq. 1281 . For cr < ctc we compare three 
possible trajectories which differ in the field at which the tran- 
sition occurs (see text for details) . The dotted line is the result 
of our numerical simulations on random graphs. 



the calculations and computed the different energy terms 
in the Hamiltonian which account for the spin-spin ex- 
change energy (?7e), the spin-random field coupling term 
{Ud), and the energy associated with the external driv- 
ing field {—HM). The analysis of the Bcthe lattice with 
coordination numbers z > 3 allows to understand (with 
analytic equations) the role played by each energy term 
in the disorder induced phase transition that separates 
the phase with smooth hysteresys loop from the phase 
with discontinuous hysteresis loop. The availability of 
the separate energy terms allows the study of the energy 
dissipation as a function of the external field along the 
hysteresis loop. 




FIG. 9: Comparison of the energy dissipated at the transition 
jump Qt (•) with the total energy dissipated along the full 
path Q_oo-.oo (■). 



In Fig. El we compare the energy dissipation associated 
with the magnetization transition jump at H2 {Qt — 
AH) with the total energy dissipated along the full path 
Q_oo^oo- One can see that for a < ac the dissipation at 
the transition Qt represents a large fraction of the total 
dissipation. 

VI. SUMMARY AND CONCLUSIONS 

We have studied the RFIM on a Bethe lattice with a 1- 
spin-flip local relaxation metastable dynamics, following 
the method proposed in Refs. 0, 0. We have extended 
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